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NUMERICAL  SIMULATIONS  FO  THE  STRUCTURE  OF 
SUPERSONIC  SHEAR  LAYERS 


1.  INTRODUCTION 

Supersonic  shear  layers  are  inherently  more  stable  than  their  subsonic  counterparts,  and 
this  results  in  the  practical  consequence  that  there  is  less  mixing.  Adequate  mixing  in  supersonic 
flows  is,  however,  an  essential  requirement  for  propulsion  devices  designed  to  provide  chemical 
reactions,  and  heat  release  in  supersonic  streams  of  air  and  fuel  For  this  reason,  an  understanding 
of  the  behavior  of  supersonic  shear  layers  and  jets  is  a  necessary  prerequisite  to  setting  criteria  for 
the  control  of  mixing  and  energy  release  in  supersonic  flows. 

This  paper  presents  a  numerical  investigation  of  the  role  of  selected  inlet  flow  variables 
on  mixing  and  stability  in  a  bounded,  unforced  supersonic  shear  layer.  More  specifically, 
two-dimensional  time-dependent  simulations  of  the  mixing  of  two  supersonic  parallel  streams  of 
gas  (initially  separated  by  a  splitter  plate)  are  performed  to  determine  the  effects  of  the  density, 
velocity,  and  pressure  differences  between  the  streams.  The  paper  reports  and  compares  the  relative 
qualitative  differences  in  mixing  diat  occurs  and  the  types  of  coherrait  structures  that  are  formed. 

In  these  simulations,  we  have  taken  particular  advantage  of  the  fact  that  computations  can 
extend  the  work  of  analysis  into  the  nonlinear  regime  for  even  broader  ranges  of  parameters,  and 
yet  the  flow  additions  are  still  idealized  and  easily  specified  and  varied.  Riysical  properties  can  be 
measured  and  shocks  and  other  flow  structures  can  be  detected,  all  of  which  are  difficult  to  observe 
experimentally.  As  with  idealized  analytical  solutions  and  basic  fluid  experiments,  the  purpose  is  to 
gain  insight  into  some  of  the  fundamental  parameters  controlling  the  behavior  of  the  flow  and 
present  guidelines  for  when  to  expect  more  mixing  and  for  choice  of  sensible  but  very  expensive 
three-dimensional  computations. 

Manuicripl  approved  April  17,  1990. 
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n.  BACKGROUND 


Subsonic  shear  layers  have  attracted  investigators  for  many  years.  Brown  and  Roshko  ^ 
made  a  fundamental  breakthrough  by  discovering  that  shear  layers  form  discrete  large-scale 
structures.  These  structures  can  become  relatively  large  and  roll  up  in  a  coherent  manner  that  greatly 
increases  the  mixing  surface.  It  was  later  shown  that  the  growth  rate  of  the  layer  is  governed  by  the 
pairing  of  these  structures.  Subsonic  jets  and  shear  layers  are  naturally  unstable  and  usually  lead  to 
this  large-scale  mixing.  The  higher  the  Mach  number,  the  longer  it  takes  for  a  jet  or  a  shear  layer  to 
become  unstable,  thus  reducing  the  amount  of  mixing  in  a  given  time  and  distance.  At  high  Mach 
numbers,  the  mixing  is  substantially  reduced.  With  growing  interest  in  supersonic  propulsion 
systems  and  hypersonic  airbreathing  vehicles,  there  is  an  immediate  need  to  understand  the  basic 
mechanisms  causing  the  reduced  mixing  and  for  overcoming  the  apparent  limits.  The  physics  of 
supersonic  shear  layers  is  presently  under  intense  research^"^. 


In  subsonic  shear  layers,  the  basic  method  of  relating  the  different  properties  of  the 
streams  is  to  normalize  with  respect  to  the  high  velocity  stream.  In  supersonic  mixing,  however, 
large  changes  in  the  density  ratio  can  cause  the  low  Mach  number  stream  to  have  a  greater  velocity 
than  the  high  Mach  number  stream.  The  basis  for  normalization  is  thus  not  clearly  defined  for 
supersonic  shear  layers.  Recent  studies^'^  have  attempted  to  use  the  'convective  Mach  number'  as  a 
normalization  parameter  in  one  way  or  the  other.  The  definitions  of  the  upper  and  lower  convective 
Mach  numbers  (for  equal  pressure  cases)  are^ 


U,  -  Uc 

M.,.  = 


U.  -  u. 


(1) 


where  Ui  and  U2  are  the  free  stream  velocities  of  the  upper  and  lower  streams  and  is  the 
convective  velocity  of  the  structures  (assumed  to  be  constant) .  The  quantities  a^  and  a2  are  the 
sonic  velocities  of  the  two  streams  at  the  inlet  condition.  For  the  overexpanded  (faster  stream  at  a 
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lower  inessure)  and  underexpanded  cases,  the  inlet  parameters  are  replaced  by  the  conditions  of  the 
two  streams  where  they  reach  static  pressure  equilibrium  after  undergoing  rarefaction  and 
compression . 


The  relation  between  and  is  traditionally  obtained  by  requiring  that  the  total 
pressure  of  the  two  streams  in  the  convective  frame  be  equal^.  This  originates  from  a  well-known 
subsonic  flow  argument,  which  state  that  dioe  exists  a  stagnation  point  between  two  structures  that 
must  be  stable,  so  that  the  pressures  at  that  point  must  balance.  Equating  the  stagnation  pressures  of 
the  two  streams  produces  an  estimate  of  the  convective  velocity  of  the  structure.The  above  argument 
was  subsequently  applied  to  supersonic  flows  by  Papamoschou  and  Roshko^,  who  assumed  that 
the  flow  comes  to  rest  isentropically  at  the  stagnation  point  Hence,  no  shock-wave  losses  were 
accounted  for.  If  the  upper  and  lower  streams  both  have  the  same  ratio  of  specific  heats,  then  the 
velocity  of  the  convected  disturbances  is  given  by 


^  «lU2  * 

a,  ♦  Sj 

and  there  is  just  one  convective  Mach  numb^ , 


(2) 


Me  =  M  .  =  M 
^  C.1  ca 

Recent  measurements  of  Papamoschou^  have,  however,  shown  that  the  convective  Mach  numbers 
corresponding  to  each  side  of  a  shear  layer  (MqJ  and  M^^)  are  very  different  thus  contradicting 
the  above  isentropic  model  of  structures  which  predicts  that  they  are  equal  or  very  close. 


Spark  Schlieren  photogrtq)hs^~^  have  shown  that  the  growth  rate  of  supersonic  shear 
layers  is  significantly  less  than  that  of  incompressible  shear  layers  layers  with  the  same  velocity 
and  density  ratios  between  two  mixing  streams.  The  drastic  reduction  in  growth  rate  has  been 
correlated  to  M^  as  defined  above.  The  experimental  measurements^  show  that  the  ratio  of  the 
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growth  rate  of  a  supersonic  shear  layer  to  the  growth  rate  of  an  incompressible  shear  layer,  (with 
the  same  velocity  and  density  ratios),  decreases  as  Mg  increases.  This  growth  rate  decreases  rapidly 
in  the  range  Me  -  0.4  -  0.5.  The  relative  reduction  again  levels  off  as  Mg  exceeds  unity.  Recent 
measurements  in  supersonic  shear  layers^"^  also  confirm  this  observation.  Linear  inviscid  stability 
analysis  has  shown^  a  similar  correlation  between  the  ratio  of  maximum  growth  rates  of  stability 
waves  of  supersonic  and  incompressible  shear  layers  and  based  on  a  frame  of  reference 
traveling  with  the  instability  waves.  For  very  small  values  of  M^,  the  reported  ratio  of  the  growth 
rate  of  a  supersonic  shear  layer  to  the  growth  rate  of  an  incompressible  shear  layer  is  near  unity. 

Experimental  studies  of  the  details  of  the  mixing  layer  in  supersonic  shear  layers  are 
difficult  as  the  mixing  layer  is  complex  when  the  streams  are  supersonic^*^.  In  addition,  the 
properties  of  supersonic  shear  layers  are  sensitive  to  small  variations  in  the  density  ratio  of  the  two 
streams .  This  in  turn  can  affect  the  optical  diagnostic  systems  used  to  characterize  the  shear  layers. 
In  such  physical  systems,  numerical  computations  often  provide  missing  information  and  behavior 
trends.  Th^  have  an  advantage  in  that  there  are  no  interference  effects  present  due  to  measurement 
probes.  In  addition,  they  can  simulate  situations  where  the  flow  conditions  are  difficult  or 
expensive  to  achieve  in  the  laboratory.  However,  as  widi  any  tool,  there  are  important  limitations  to 
what  the  simulations  can  do.  The  validity  of  the  results  depends  on  the  computational  domain 
chosen,  the  particular  algorithm  used  to  solve  the  model  equations  and  its  implementation  of  the 
initial  and  boundary  conditions,  the  particular  computer,  and  even  the  interpretation  of  the  output 
data. 


Numerical  simulations  have  been  and  are  currently  a  powerful  tool  for  investigating  the 
evolution  of  subsonic  shear  layers^^^^.  Due  to  the  presence  of  complicated  shock  structures  and 
sharp  gradients,  accurate  simulations  of  the  supersonic  shear  layers  warrant  algorithms  with 
minimal  numerical  diffusion.  Recently  Guirguis  et  aL^  have  demonstrated  the  applicability  of 
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numerical  simulations  to  the  study  of  fundamental  processes  in  time-dependent  two-dimensional 
supersonic  shear  layers.  The  Flux-Corrected  Transport  (FCT)  algorithm  was  used  to  study  the 
mixing  and  structure  in  supersonic  shear  layers.  Soestrino  et  al.^^  studied  the  evolution  of 
instabilities  and  the  physical  processes  associated  with  turbulent  mixing  of  temporally  supersonic 
shear  layers,  i.e.,  shear  layers  with  periodic  boundary  conditions  in  the  flow  directions.  A 
second-order  total  variation  diminishing  (TVD)  algorithm  was  used.  Numerical  simulations  of  shear 
layers  in  open  domains  were  recently  carried  out  by  Lele^^.  In  this  case,  the  two  streams  were 
matched  in  static  pressure  but  had  different  densities.  The  spatially  evolving  simulations  presented 
were  forced  at  the  inflow  by  adding  a  small  normal  velocity  disturbance  at  the  most  unstable 
frequency  and  its  first  two  subharmonics.  In  the  temporally  evolving  simulations,  small  amplitude 
incompressible  disturbances  were  added  to  the  iititial  tangent  hyperbolic  flow  profile.  The  behavior 
of  a  plane  free  shear  layer  was  studied  by  Tang  et  al.^^.  A  second-order  ADI  procedure  as  well  as 
a  modified  McCormack  algorithm  were  used  to  obtain  the  solutions.  Small-ampUtude  oscillations  in 
the  normal  velociQr  were  found  to  grow  as  they  convected  downstream. 

This  work,  following  Guirguis  et  al.^,  is  a  first  attempt  to  do  a  systematic  numerical 
parameter  study  of  unforced  and  confined  supersonic  shear  layers.  The  procedure  is  to  select  a  base 
case  and  vary  the  parameters  around  this  case  in  a  systematic  way.  The  base  case  consists  of  two 
(equal  pressure)  parallel  streams  of  supersonic  gas  at  inlet  Mach  numbers  of  4.0  and  1.5, 
respectively,  flowing  in  a  chamber  bounded  on  the  sides.  The  density,  velocity,  and  pressure  in  the 
two  streams  are  varied  so  that  equal  pressure,  underexpanded,  and  overexpanded  systems  are 
considered.  In  the  calculations  presented  hoe,  however,  the  effective  momentum  thickness  is 
generally  held  fixed. 

The  present  work  extends  the  previous  work  in  several  ways.  First,  the  parameters  were 
systematically  varied  to  look  for  the  regimes  of  enhanced  mixing.  The  evolution  of  the  major 
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physical  quantities  in  an  unforced  system  is  presented  and  these  are  used  to  analyze  the  dynamics 
of  the  flow  and  the  mixing  characteristics.  Two  diagnostics,  flow  visualization  and  fourier 
analysis,  were  used  in  this  study  to  help  analyze  the  computed  results.The  convective  velocities 
MqJ  and  are  directly  measured  from  the  predictions  and  compared  to  the  isentropic  value. 
Since  the  same  medium  (air)  was  considered  for  both  streams,  for  equal  pressure  cases,  the  density 
variation  of  the  streams  was  caused  by  temperature  differences. 

m.  THE  NUMERICAL  MODEL 

The  physical  system  we  simulate,  shown  schematically  in  Fig.  1,  is  a  confined  mixing 
layer  formed  by  two  streams  of  gases  iiutially  separated  by  a  splitter  plate.  The  flow  domain  is 
bounded  by  two  perfectly  smooth  parallel  walls.  The  splitter  plate  is  not  included  in  the 
computational  domain  and  its  trailing  edge  forms  the  left  boundary  of  the  computational  domain. 
For  supersonic  essentially  inviscid  flow  calculations,  this  does  not  introduce  significant  errors  as 
there  is  no  feedback  of  information  upstream.  The  initially  laminar  streams  have  a  top-hat  axial 
velocity  profile  at  the  trailing  edge  of  the  splitter  plate  and  this  profile  is  held  constant  during  the 
computations.  The  calculations  are  unforced  in  the  sense  that  no  external  frequency  is  imposed  on 
the  flow.  These  are  idealized  initial  conditions  since  most  laboratory  tunnel  flows  are  not  laminar  at 
the  Mach  numbers  considered  in  this  paper . 

For  the  calculations  rqiorted,  the  inflowing  upper  stream  is  always  supersonic  and  at 
a  higher  velocity  than  the  inflowing  lower  stream.  The  lower  stream  can  be  either  supersonic  or 
substHuc .  By  analogy  to  supersonic  jets,  an  undoexpanded  shear  layer  is  defined  as  the  one  in 
which  the  pressure  of  the  faster  stream  is  higher  at  the  entrance  than  that  of  the  slower  stream  and 
vice  versa  for  overexpanded  shear  layers.  Underexpanded,  equal-pressure,  and  overexpanded  shear 
layera  were  considered.  As  shown  in  Figure  1,  properties  of  the  faster  stream  are  denoted  by  a 


6 


subscript  1,  and  the  slower  flow  by  2. 


The  numerical  model  used  to  carry  out  the  simulations  consists  of  the  two-dimensional, 
time-dependent  inviscid  conservation  equations  of  mass  density ,  momentum  and  energy  density  for 
an  ideal  gas: 


ap 

at  ” 

-V.(pV) 

(3) 

a(pv) 

at 

=  -V.(pW)  -  VP 

(4) 

dE  _ 

at 

-V.(EV)  -  V.PV 

(5) 

Here  E  =  P/  (y  -  1)  +  (1/2)  pV2  is  the  total  energy  and  V  (where  V  =  iU  +  jV),  P,  p  and  y 
axe  the  fluid  velocity  vector,  pressure,  total  mass  density,  and  the  ratio  of  the  specific  heats, 
respectively.  To  compute  the  mixing  (as  a  result  of  convection),  each  of  the  two  streams  is 
convected  separately,  allowing  us  to  calculate  the  mixing  in  any  given  computational  cell.  This  was 
done  by  considering  the  conservation  equations  for  individual  species  densities  for  the  upper  and 
lower  streams. 

The  equations  are  solved  using  the  nonlinear,  high-order,  explicit,  compressible  finite- 
difference  FCT  algorithm  (Boris  and  Book^^).  Through  a  two-step  predictor  corrector  scheme, 
FCT  ensures  that  all  conserved  quantities  remain  monotonic  and  positive.  First,  it  modifies  the 
linear  properties  of  a  high  order  algorithm  by  adding  diffusion  during  convective  transport  to 
prevent  dispersive  ripples  from  arising.  The  added  di^sion  is  then  removed  in  an  antidiffusion 
phase  of  the  integration  cycle.  Hence  these  calculations  maintain  the  high  order  of  accuracy  without . 
requiring  artificial  viscosity  to  stabilize  the  algorithm.  The  algorithm  has  been  extensively  tested 
and  used  in  the  last  ten  years  (see,  for  example  bibliogr!q>hy  in  Oran  and  Boris^^)  to  predict  a  wide 
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variety  of  reactive  and  nonreactive  flows.  The  explicit  FCT  algorithm  was  used  (Oran  et  al.^8)  to 
investigate  the  properties  of  condensational  instabilities  in  solar  and  astrophysical  plasmas,  (jiowth 
rates,  decay  rates  and  oscillation  frequencies  of  the  perturbations  determined  from  the  linear  analysis 
were  in  excellent  agreement  with  the  simulations.  In  addition,  the  FdT  algorithm  has  been  used  for 
the  analysis  of  a  variety  of  plasma  and  fluid  dynamic  instability  problems^  Employing  the 

FCT  algorithm,  Grinstein  et  al.^^~^^  calculated  the  vortex  dynamics  and  asymmetric  mixing  in 
compressible  subsonic  plane  and  axisymmetric  shear  layers  in  both  two  and  three  dimensions. 
Also  using  this  algorithm,  Kailasanath  et  al.^^  have  calculated  the  behavior  of  cold  and  chemically 
reacting  axisymmetric  confined  jets.  The  present  calculations  were  carried  out  using  the  code 
LCPFCT  which  uses  the  most  recent  one-dimensional  version  of  Flux-Corrected  Transport 
technique  with  fourth-order  accurate  phases  and  minimum  residual  diffusion.  A  timestep-splitting 
technique  is  used  to  couple  the  two  dimensions.  The  calculations  are  vectorized  to  take  full 
advantage  of  the  architecture  of  the  CRAY-XMP  supercomputer  used  for  the  calculations. 

No  explicit  subgrid  turbulence  model  was  included.  However,  there  is  an  effective  filter 
for  high  frequencies  which  is  part  of  the  nonlinearity  of  the  FCT  algorithm.  The  filter  does  not 
affect  die  large-scale  structures,  but  ensures  that  wavelengths  smaller  than  a  few  computational  cells 
are  numerically  diffused.  Thus,  provided  that  the  computational  grid  is  fine  enough  to  resolve  the 
large-scale  features  of  the  flow,  the  residual  numerical  viscosity  of  the  algorithm  mimics  the 
behavior  of  physical  viscosity  at  high  Reynolds  number  by  damping  small-scale  structures  on  the 
order  of  a  few  computational  cells.  There  are  essentiaUy  no  viscous  effects  at  scales  greater  than 
four  or  five  computational  cells. 

A25cmx3cm(XxY)  computational  domain  was  considered  for  the  majority  of  the 
cases  studied.  The  splitter  plate  was  located  symmetrically  between  the  confining  walls.  The  two 
streams  are  ideal  gases  (air)  with  constant  specific  heats.  Inflow  and  outflow  boundary  conditions 
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must  be  used  to  represent  the  effects  of  the  outside  world  on  the  computational  domain.  Ideally, 
such  boundary  conditions  would  not  allow  unphysical  reflections  or  wave  absorption  at  the 
boundaries  of  the  computational  domain.  Imposing  such  boundary  conditions  is  difficult 
numerically  and  is  currently  the  topic  of  a  great  deal  of  numerical  analysis  research.This  difficulty  is 
most  evident  in  numerical  simulations  of  subsonic  flows,  but  it  is  not  such  a  major  problem  in 
supersonic  flows.  The  inlet  conditions  of  the  dependent  variables  were  specified.  Zero-gradient 
outflow  conditions  (at  X  =  L)  were  used  for  all  variables.  The  top  and  bottom  of  the  computational 
region  (Y  =  0  and  H)  are  perfectly  reflecting  hard  walls.  A  191  x  61  (X  x  Y )  imiform  grid  was 
used  for  most  of  the  calculations  reported.  The  mesh  was  uniform  in  order  to  ensure  that 
downstream  and  deflected  structures  were  well  resolved.  A  series  of  calculations  were  performed 
where  the  inlet  Mach  numbers,  velocity ,  pressures  and  densities  of  the  two  streams  at  the  inlet 
were  systematically  varied.  Extended  domains  with  the  same  grid  density  ,40  cm  x  3  cm  and  40 
cm  X  5  cm  and  denser-grid  calculations  (382  x  122  and  382  x  61)  on  the  same  25  cm  x  3  cm 
domain  were  also  performed  for  a  specific  base  case.  The  Courant  number  was  held  at  0.5  to 
insure  numerical  stability  and  to  resolve  the  time  evolution  properly.  This  condition  was  based  on  a 
velocity  equal  to  sound  speed  plus  flow  speed  characteristic  of  the  faster  stream. 

The  boundary  layers  and  the  actual  velocity  profiles  on  the  splitter  plate  are  important 
parameters  since  they  provide  a  characteristic  transverse  length  scale  for  the  shear-layer  instability. 
Since  different  experimental  set-ups  have  different  boundary  layers  and  velocity  profiles  developing 
on  the  splitter  plate  the  characteristic  length  scale  Sq  (the  initial  momentum  thickness)  can  vary 
from  one  set  of  experiments  to  another.  Our  calculations  were  started  with  a  step-function  velocity 
profile  (with  a  slope  AU/AY  where  AU  =  Uj  -  U2  and  AY  is  the  cell  width  in  the  Y  direction) 
across  the  shear  layer,  i.e.,  an  infinitesimally  thin  vortex  sheet  spanning  the  whole  domain,  from 
the  trailing  edge  of  the  splitter  plate  up  to  the  right  boundary.  These  profiles  smooth  out  during  the 
initial  timestqps  and  relax  to  a  shear-layer  profile  with  an  effective  initial  momentum  thickness  Qq. 
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As  found  in  experiments,  Qq  can  vary  among  diffo’ent  simulations  performed  for  example,  using 
different  algorithms.  Therefore  it  is  important  to  characterize  the  relevant  lengthscales  in  either 
calculations  or  experiments.  Earlier  calculations  using  the  approach  adopted  in  this  paper  have 
shown  that  a  step-profile  for  the  inlet  velocity  produces  an  estimate  of  Qq  as  one  or  two  cell-width 
in  the  cross  stream  direction^.  Our  present  calculations  also  support  this  observation,  which  means 
that  we  can  detrarnine  an  effective  initial  momoitum  thickness  by  setting  the  grid  spacing. 

The  small  perturbation  which  initiates  the  transition  from  laminar  to  transitional  flow 
occurs  at  the  first  time  step  of  the  calculations.  The  perturbation  generates  small  pressure  gradients 
and  diffuses  vorticity  at  the  shear  layer,  near  the  edge  of  the  splitter  plate.  Due  to  the  smoothing  of 
the  initial  step-function  profrle,  the  initially  uniform  pressure  field  (for  equal  pressure  cases) 
develops  a  small  gradient  across  the  thickening  shear  layer  in  the  region  close  to  the  edge  of  the 
splitter  plate  where  the  free  streams  meet  This  disturbance  moves  downstream  as  the  integration 
proceeds,  g»ierating  the  transverse  flows  which  trigger  the  instability. 

The  numerical  simulations  predict  values  of  the  density,  momentum,  energy  and 
number-density  ratio  for  each  of  the  computational  cells  as  a  function  of  time.  The  evolution  of  the 
flow  is  shown  through  sequences  of  contours  of  number-density  ratio,  where 

R  =  Ni/(Ni  +  N2) 

TheNj  and  N2  are  the  number  densities  corresponding  to  the  top  and  bottom  inflowing  streams. 
The  value  of  R  varies  between  0  and  1;  a  value  of  0  cOTiesponds  to  100%  of  the  bottom  stream  and 
a  value  of  1  corresponds  to  100%  of  the  top  stream. 

Kailasanath  et  al.^^  have  shown  that  the  peaks  in  velocity-  and  pressure-fluctuations 
spectra  taken  at  fixed  locations  can  be  related  to  the  organization  and  coherence  in  the  flow.  The 
fluctuations  at  a  given  location  reflect  the  local  passage  of  large  structures  as  well  as  information 
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about  events  at  other  locations  of  the  flow  transmitted  by  sound  waves,  such  as  vortex  mergings 
that  generate  fluctuations  in  compressible  flows.  A  strong  narrow  spectral  line  usually  indicates  that 
a  structure  passes  the  location  at  a  very  regular  interval.  We  have  therefore  performed  fourier 
analysis  of  pressure  and  velocity  fluctuations  at  selected  positions  in  the  computational  domain  to 
help  determine  the  organization  and  coherence  in  the  flow. 

IV.  RESULTS  and  DISCUSSION 

The  computations  performed  are  divided  into  three  broad  categories: 

1.  Both  streams  are  at  the  same  pressure  and  are  supersonic  as  they  pass  the  tip  of  the  splitter  plate. 
Calculations  were  performed  by  either  keeping  the  density  uniform  and  varying  the  velocity  ratios 
or  by  maintaining  the  velocity  ratio  constant  and  varying  the  density  ratios. 

2.  Both  streams  are  at  the  same  temperature  and  supersonic  at  the  inlet,  but  the  pressures  of  the  two 
streams  at  the  inlet  are  different  Both  underexpanded  and  overexpanded  cases  were  considered. 

3.  Both  streams  are  at  the  same  pressure  and  temperature,  but  the  upper  stream  is  supersonic  and 
the  lower  stream  is  subsonic. 

The  convective  Mach  number,  ,  as  obtained  by  equation  (2)  was  computed  for  each  case  and  are 
reported  with  the  inlet  parameters. 

1.1  Equal-Pressure  Streams 
a.  Uniform  density 

The  convective  Mach  number  and  the  velocity  ratios  were  systematically  varied  by 
decreasing  the  Mach  number  of  the  lower  stream ,  M2 ,  while  keeping  the  Mach  number  of  the 
upper  stream,  Mj,  constant  Conditions  for  the  base  case  are  given  in  Table  I.  For  these 
conditions,  the  value  of  Mq  is  directly  proportional  to  (M^  >  M2),  so  that  M^  =  125  . 
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Figure  2  shows  the  contours  of  the  mixing  ratio  R  within  the  computational  domain  at  a 
sequence  of  times  corresponding  to  8.6  ps  intervals  or  SO  time  steps,  at  a  range  of  time  well  past 
the  initial  stages  of  the  development  of  the  flow  in  the  calculations.  The  shear  layer  initially 
becomes  unstable  at  approximately  9  cm  from  the  tip  of  the  splitter  plate.  Subsequently  coherent 
structures  develop  and  move  downstream.  The  paths  of  some  of  these  organized  structures  in  the 
figure  show  that  there  is  merging  indicating  that  the  convective  velocities  of  all  the  structures  are 
not  the  same.  The  convective  velocities  of  these  structures  (as  computed  from  Figure  2) 
significantly  differ  from  the  (theoretical)  value  of  « 125.  For  instance,  fee  convective  velocity 
of  a  typical  structure,  the  one  at  fee  left  undergoing  merging  wife  a  smaller  structure,  was 
11.6  X  10^  cm/s  which  gives  -  0.67  and  Mc^  =  1.83.  As  expected,  this  does  not  agree 
wife  predictions  of  fee  fee  isentropic  model  The  above  results  are,  however,  consistent  wife  fee 
recent  experimental  observations  by  Papamoschou^. 

Figure  3  shows  the  density  contours  of  the  flow  field  at  fee  end  of  time  step 
4000.  The  density  field  remains  uniform  until  the  shear  layer  becomes  unstable.  A  system  of 
oblique  shock  waves  is  observed  attached  to  fee  large  structures  in  fee  lower  half  of  fee  domain 
where  fee  computed  =  1.83.  The  evolution  of  the  flowfield  shown  in  Figure  2  can  be  related 
to  the  fouiier  analysis  of  fee  pressure  and  velocity  fluctuations  observed  at  specific  downstream 
locations  in  the  computational  domain.  In  Figure  4  we  show  fee  spectra  of  fee  transverse  velocity 
and  pressure  fluctuations  at  the  point  (18.75  cm,  2.0  cm),  indicated  by  an  "X"  in  Figure  1.  The 
location  was  chosen  after  sanq)le  calculations  were  made  and  regions  wife  distinctive  structures 
were  identified.  Both  fee  pressure  and  velocity  fluctuations  show  a  dominant  frequency  of  30  kHz, 
corresponding  to  fee  passage  of  the  structures  shown  in  Figure  2.  The  fourier  analyses  do  not 
show  strong  specific  subharmonics  signifying  that  fee  passing  structures  do  not  merge  in  any 
clearly  synchronized  way. 
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Table  I.  Base  case  conditions 


1 

2 

p 

1.17  X  10*3  gm/cm^ 

1.17  X  10"3  gm/cm^ 

T 

300K 

300K 

P 

1.01  X  lO^dynes/cm^ 

1.01  X  10^  dynes/cm3 

U 

13.95  X  10^  cm/s 

522  xlO^  cm/s 

M 

4.0 

15 

Figure  5  compares  contours  of  the  mixing  ratio  R  at  timestep  4000  for  calculations  in 
which  M2  was  gradually  increased,  M2  -  2.0,  3.0,  and  3.S.  The  value  of  and  all  other 
conditions  were  kept  the  same  as  the  base  case,  which  resulted  in  =  1.0, 0.5,  and  0.25  and 
velocity  ratios  U1AJ2 -2.0, 1.33,  and  1.13,  respectively.  The  contours  ofR  show  that  as  the 
velocity  ratio  decreases,  the  shear  layer  spreads  less  and  becomes  unstable  further  downstream,  a 
result  also  seen  in  subsonic  shear  layers .  The  spectra  of  the  pressure  and  velocity  fluctuations 
become  noisier  as  and  the  velocity  ratio  decrease  widi  a  proliferation  of  peaks  between  40  kHz 
and  80 1^  range  in  the  spectrum. 

Calculations  were  also  carried  out  widi  Mj  decreased  to  2.5  and  M2  remaining  the  same 
as  in  the  base  case.  This  corresponds  to  a  velocity  ratio  U1/U2  1.66  and  M^^  %  0.5.  The  contours 
of  R,  shown  in  Figure  6  at  the  end  of  timestep  4000  ,  show  a  very  stable  shear  layer  with  no 
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convective  mixing.  The  slight  spread  in  the  shear  layer  is  caused  by  the  small  though  nonzero 
residual  numerical  diffusion  associated  with  the  finite^ifference  scheme.  Figures  5(b)  and  S(c) 
show  that  the  shear  layer  becomes  unstable  in  the  same  computational  domain  for  cases  where  the 
velocity  ratio  U  i  AJ2  is  lower  than  1.66.  In  particular,  the  value  of  is  same  for  cases  shown  in 
Figures  5(b)  and  6.  However,  Mi  =  4.0  and  M2  =  3.0  in  Figure  5(b)  while  the  corresponding 
values  in  Figure  6  are  2.5  and  1.5  respectively.  Since  the  calculations  were  carried  out  without  any 
forcing,  the  pressure  fluctuations  are  unable  to  trigger  the  Kelvin-Helmholtz  instability  within  the 
length  of  the  computational  domain  in  Figure  6.  Guirguis  et  al.^  have  shown  earlier  that  an 
unforced  supersonic  shear  layer,  that  appears  stable  in  a  short  computational  domain,  may 
eventually  become  unstable  further  downstream. 

To  study  the  flow  in  a  longer  domain,  the  following  cases  were  considered  with  the  inlet 
conditions  identical  to  the  base  case.  First,  the  computations  were  carried  out  on  a  domain  with 
dimensions  of  40  cm  x  3  cm  with  a  300  x  61  grid.  This  gave  the  same  grid  density  as  in  the  base 
case  but  the  domain  is  15  cm  longer.  Figure  7  depict  the  contours  of  R  for  this  case  at  the  end  of 
time  step  4000.  Here  the  growth  of  the  naturally  unstable  shear  layer  was  obstructed  downstream 
by  the  walls.  The  domain  dimensions  were  then  increased  to  40  cm  x  5  cm  with  a  300  x  101  grid 
to  keep  the  same  grid  density.  Figure  8  shows  contours  of  the  mixing  ratio  R  within  the 
computational  domain  at  a  sequence  of  times  for  the  extended  domain.  As  in  Figure  2,  the  frames 
in  Figure  8  are  8.6  ps  apart  which  corresponds  to  50  timestep  intervals.  As  in  the  base  case,  the 
instability  is  triggered  near  9  cm  .  However,  the  mixing  increases  substantially  in  the  extended 
length  of  the  domain.  The  slight  differences  in  the  contours  in  the  first  25  cm  of  Figure  8  and 
Figure  2  are  attributed  to  the  different  heights  of  the  confining  chambers,  which  in  turn  impose 
slightly  different  acoustic  properties  on  the  system.  The  computed  velocities  for  M^^i  and  M^  2 
(for  a  typical  structure)  were  close  to  those  found  from  Figure  2.  These  calculations  also  show  that 
in  this  study,  the  location  of  the  outflow  boundary  has  little  effect  on  the  flow  field  predictions 
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upstream. 


For  simulations  performed  by  nonlinear  monotone  methods  such  as  FCT,  the  residual 
numerical  diffusion  is  finite  but  very  low.  Hence  for  simulations  that  do  not  include  an  explicit 
model  for  physical  viscosity  or  for  which  die  jidiysical  viscosity  is  less  than  the  numerical  viscosity, 
the  elective  Reynolds  numbers  for  various  sized  structures  in  the  flow  are  determined  by  the  mesh 
size.  Because  of  this  inherent  nonlinearity  of  the  method,  the  effective  Reynolds  number  varies 
according  to  the  size  of  the  structure  resolved  relative  to  the  mesh  spacing.  For  example  if  a 
structure  is  resolved  by  twenty  computational  cells,  the  effective  Reynolds  number  is  very  high  and 
the  flow  is  essentially  inviscid  for  this  structure.  If  the  structure  is  resolved  by  three  cells,  the 
effective  Reynolds  number  is  low  and  the  structure  can  be  dissipated.  The  nonlinearity  of  the 
method  means  that  the  progression  to  increased  effective  Reynolds  number  is  faster  than  linear. 
Thus  when  mesh  size  is  decreased,  smaU-sxale  structures  of  the  flow  are  resolved  as  long  as  the 
physical  viscosity  is  less  than  the  numerical  viscosity.  For  larger  mesh  size,  the  small  scales  are 
filtered  out  In  order  to  study  the  effects  of  grid  size  and  the  inlet  velocity  profile  on  the  solutions, 
the  base  case  calculations  woe  done  with  different  computational  grids,  viz.  382  x  61  (twice  as  fine 
in  the  X  direction)  and  382  x  122  (twice  as  fine  in  both  X  and  Y  directions).  Contours  of  mixing 
ratio  R  at  the  end  of  timestep  KXX)  is  shown  in  Figures  9(a)  and  9(b)  for  the  base  case  ( 191  x  61) 
and  for  the  382  x  61  grid  respectively.  For  these  two  cases,  the  inlet  velocity  profile  ( and  hence 
)  is  idratical  as  the  cross  stream  cell  widfo  is  the  same.  Though  the  essential  large  scale  features 
shown  in  the  two  figures  are  similar,  more  small-scale  structures  are  found  in  Figure  9(b)  due  to  the 
finer  grid  size  in  the  X  direction  only.  Figure  9(c)  shows  the  mixing  contours  obtained  with  a  grid 
of  382  X  122  at  the  same  physical  time  and  precisely  the  same  inlet  velocity  profile  considered  in 
Figures  9(a)  and  9(b).  Due  to  the  doubling  of  die  grid  in  the  Y  direction,  AY  is  half  the  size  used  in 
the  above  cases.  The  inlet  velocity  profile  was,  however,  spread  linearly  across  two  cross  stream 
cells  so  that  it  was  identical  to  those  applied  for  Figures  9(a)  and  9(b).  The  agreement  between  the 
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predictions  shown  in  Figures  9(a)  and  9(c)  is  good.  As  the  solution  progressed  in  time,  the 
similarities  of  the  above  predictions  became  less  evident  due  to  the  effective  filtering  of  high 
frequencies.  The  filter  does  not  affect  the  large-scale  structures,  but  ensures  that  wavelengths 
smaller  than  a  few  computational  cells  are  numerically  diffused.  It  is,  however,  evident  that  the 
shape  of  the  inlet  profile  plays  a  key  role  in  the  development  of  the  instability  and  for  the  present 
calculations,  the  large-scale  structures  are  essentially  grid  indq)endent 

Predictions  were  also  obtained  with  the  382  x  122  grid  where  the  inlet  velocity  profile 
was  across  one  cross  stream  cell  only.  Figures  10(a)  and  10(b)  compare  the  contours  of  R  with 
the  denser  grid  (with  the  above  inlet  velocity  profile)  to  those  of  R  for  the  base  case  at  the  same 
physical  time .  Note  that  because  Oq'  is  halved,  the  physics  of  the  two  flows  is  different.  In  the 
denser-grid  case,  the  mixing  instability  occurs  at  a  shorter  downstream  distance,  about  4  J  cm 
from  the  trailing  edge  of  the  flitter  plate,  compared  to  the  base-case  prediction  where  the  instability 
is  triggered  near  9  cm  .  However,  there  is  a  physical  similarity  in  the  two  solutions:  the  flow 
structures  in  the  ffrst  half  of  Figure  KKa)  look  very  much  like  those  in  Figure  10(b),  only  half  the 
size.  The  shedding  frequencies  computed  by  the  fotirier  analyses  of  the  pressure  and  velocity 
fluctuations  (as  shown  in  Figure  11)  were  approximately  twice  the  corresponding  shedding 
frequencies  predicted  in  the  base  case.  The  similarity  in  the  solution  of  equations  is  consistent  with 
the  fact  that  Qq  was  reduced  by  reducing  the  cross  stream  cell  width . 


b.  Constant  Velocity  Ratios 

For  the  equal-pressure  cases,  the  density  ratio  was  varied  by  changing  the  temperature 
ratios  of  the  two  streams  but  keeping  the  velocity  ratios  constant  The  velocities  of  the  upper  and 
lower  streams  were  Ui  ®  8.71  x  10^  cm/s,  U2  =  5.22  x  10^  cm/s  respectively  ,  so  that  U 1/U2 
>  1.66.  Results  were  obtained  for  two  situations,  one  in  which  the  faster  stream  is  four  times 
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denser  than  the  slower  stream  (pi  /p2  »  4.0  and  T1/T2  «  1/4 )  and  the  other  in  which  the  faster 
stream  is  four  times  lighter  than  the  slower  stream  (pi /p2  =  1/4  and  T1/T2  =4.0).  For  the 
case  where  the  faster  fluid  is  denser,  Mj  =  S.O  and  M2  =  1.5,  and  for  the  case  where  the  faster 
fluid  is  lighter,  Mj  =  2.5  and  M2  -  3.0.  It  is  interesting  to  note  that  for  both  of  these  cases,  M^  = 
0.67. 


The  instantaneous  contours  of  R  for  the  two  cases  described  above  are  shown  in 
Figures  12(a)  and  12(b)  at  the  end  of  time  step  4(X)0.  Due  to  the  temperature  changes  in  the  streams 
due  to  density  differences,  the  time  steps  must  be  slightly  different  to  satisfy  the  imposed  Courant 
number  condition  for  stability  and  accuracy.  The  case  where  the  faster  fluid  is  lighter  (hotter),  is 
found  to  be  stable  and  no  convective  mixing  occurs.  However,  significant  mixing  occurs  when  the 
faster  fluid  is  denser  (colder)  and  here  the  contours  of  R  show  organized  structures. 

It  has  been  shown  from  linear  themy  that  the  effect  of  temperature  differential  is  different 
for  subsonic  disturbances  as  opposed  to  superscmic  disturbances  (Gropengiesser^4^  Ragab  and 
Wu^).  Brown  and  Roshko^  investigated  the  density  effects  in  subsonic  turbulent  mixing  layers. 
According  to  this  study,  when  the  faster  stream  is  denser  (U1/U2  =  /7  and  Pi/P2  =  7),  decreased 
spreading  of  the  shear  layer  was  observed  compared  to  the  case  where  the  faster  stream  is  slightly 
lighter  than  the  slower  stream  (U1/U2  =  ^7  and  Pi/P2  =  0.965).  Our  present  results  for  the 
supersonic  shear  layer  and  the  observations  by  Brown  and  Roshko  for  the  subsonic  shear  layer  are 
thus  consistent  and  agree  with  the  predictions  by  the  linear  theory. 


2.11 


U'H  1-^.1 


indcd  and  OveriBTpgnded  rases 


When  the  pressures  of  the  two  streams  are  different,  the  high-pressure  stream  expands 
dirough  a  rarefaction  fan  centered  at  the  tip  of  the  splitter  plate  while  the  low-pressure  stream  is 
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compressed  to  the  same  pressure  through  an  oblique  shock  wave  also  attached  to  the  tip  of  the 
splitter  plate.  The  actual  flow  parameters  affecting  the  properties  of  the  shear  layer  are  the  values 
ahead  of  these  waves.  The  conditions  ahead  of  the  shock  or  rarefaction  waves  that  form  for  the 
overexpanded  and  underexpanded  cases  will  be  denoted  by  1'  and  2*  for  the  top  and  the  bottom 
streams,  respectively,  as  was  done  in  Guirguis  et  al.^. 

Figure  13  shows  the  density  contours  in  the  computational  domain  at  a  time  step  4000  for 
an  underexpanded  case  where  Pi  »  2.02  x  10^  dyn^cm^,  P2  «  1.01  x  10^  dynes/cm^  and  pi 
=  2.34  X  lO"^  gm/cm^,  P2  =  1.17  x  10*3  gm/cm^.  The  remaining  parameters  are  identical  to 
those  shown  in  Table  1.  For  this  case,  the  temperatures  of  the  two  stream  at  the  tip  of  the  splitter 
plate  are  equal  (300  K).  As  a  result  of  the  pressure  difference,  the  shear  layer  turns  4.5^  towards 
the  lower  wall.  The  higher  pressure  flow  on  top  expands  through  a  rarefaction  fan  centered  at  the  tip 
of  the  splitter  plate  to  a  pressure  Pj*  =  1.232  x  10^  dynes/cm^,  while  the  lower  pressure,  slower 
flow  compresses  through  an  oblique  shock  wave  to  the  same  pressure  P2'  =  1.232  x  10^ 
dynes/cm2.  Ahead  of  the  rarefaction  fan,  Mj’  =  4.38,  Tj’  =  2632  K  and  pi’  =  1.676  x  10"3 
gm/cm^.  Ahead  of  the  shock  wave,  M2'  =  1.373,  T2'  =  308  K  and  P2'  =  1.32  x  10"^  gm/cm^. 
The  (isentropic)  convective  Mach  number  at  the  above  conditions  is  M^.  s  1.36.  The  density 
cmitours  illustrate  the  interaction  between  the  reflected  rarefaction  and  ^ock  waves  at  the  walls  and 
the  shear  layer.  The  oblique  shock  wave  starting  at  the  splitter  plate  reflects  from  the  bottom  wall. 
When  the  shock  intersects  the  shear  layer,  it  is  partially  transmitted  and  partially  reflected  due  to  the 
higher  density  of  the  top  layer.  This  produces  a  perturbation  that  triggers  the  Kelvin-Helmholtz 
instability.  Further  complex  interactions  of  the  walls  with  transmitted  and  reflected  shocks  occur 
in  the  downstream  flow. 

Figure  14  shows  the  contours  of  R  at  a  sequence  of  times.  As  in  Figure  2,  each  frame 
in  Figure  14  is  8.6  ps  apart  which  corresponds  to  intervals  of  SO  time  steps.  The  time  sequences  of 
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R  exhibit  enhar''ed  mixing  and  organization  of  the  structures  compared  to  those  shown  in  Figure  2. 
This  is  perhrqrs  due  to  the  presence  of  shocks  in  the  vicinity  of  the  splitter  plate  which  perturbs  the 
shear  layer.  The  process  of  merging  is  indicated  by  the  two  dashed  lines  coming  together  at  step 
4250.  The  computed  values  of  and  for  the  structure  at  the  left  undergoing  pairing  in 
Figure  13  are  0.918  and  2.22  respectively  as  opposed  to  the  theoretical  value  of  =  1.36. 
Since  Me, 2  ^  1-0,  oblique  shocks  waves  are  thus  found  in  the  lower  half  of  the  computational 
domain  in  Figure  13. 

In  Figure  IS  we  show  the  spectra  of  the  transverse  pressure  fluctuations  at  18.75  cm 
from  the  left  boundary  and  2.0  cm  from  the  lower  wall.  Comparing  the  results  to  those  in  Figure 
4(b)  show  that  a  second  strong  peak  appears  in  the  pressure  fluctuation  spectrum  reflecting  the 
more  organized  merging  process  seen  in  Figure  14. 

In  order  to  investigate  further  the  effects  of  underexpansion  and  overexpansion  on  the 
mixing  behavior,  the  inlet  parameters  for  die  case  shown  in  Figure  6  (that  exhibited  little  convective 
mixing)  was  considered.  Both  underexpanded  (P1/P2  =  2.0,  Pi/P2  =  2.0)  and  overexpanded  cases 
(P1/P2  =  0.5,  Pi/P2  =  0.5)  were  considered.  For  both  cases,  Mj  =  2.5  and  M2  =  1.5  and  Tj  = 
T2  =  300  K.  Figures  16(a)  and  16(b)  show  the  contours  of  R  at  a  given  time  (timestep  40(X))  for 
the  underexpanded  and  the  overexpanded  cases.  For  the  underexpanded  case,  the  shear-layer 
deflection  is  6.0^  while  in  the  overexpanded  cases,  the  deflection  is  about  5°.  These  cases  show 
that  both  underexpansion  and  overexpansion  enhance  mixing  compared  to  the  equal-pressure  case 
shown  in  Figure  6.  However,  convective  mixing  in  the  overexpianded  case  is  significantly  less  than 
the  underexpanded  case  where  the  inlet  temperatures  of  the  two  streams  are  the  same.  Both  cases 
show  organized  large-scale  mixing  and  the  structures  appear  similar  in  shape  to  those  observed  in 
subsonic  shear  layers.  The  shear  layer  becomes  unstable  near  15  cm  from  the  trailing  edge  (for  the 
underexpanded  case)  of  the  splitter  plate.  The  velocity  and  density  ratios  in  the  shear  layer 
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immediately  ahead  of  the  shock  and  the  rarefaction  waves  are  V\f\52  ^  1*96  and  p^V  P2  = 
126  for  the  underexpanded  case  and  Ui'/U2'  =  1.45  and  pi'/  P2'  =  0.81  for  the  overexpanded 
case. 


3.  SupcggnigiSub»nig  Ii 


iS 


A  numbo*  of  experimental  studies  have  been  performed  where  a  supersonic  flow  meets 
subsonic  flow  at  the  trailing  edge  of  a  splitter  plate  4^  jg  section,  the  results  of  a  limited 
computational  study  carried  out  on  confined  shear  layers  formed  by  a  supersonic  and  a  subsonic 
stream  are  reported.  The  mathematical  formulation  and  the  boundary  conditions  were  identical  to 
those  reported  in  an  earlier  section.  The  computadonal  domain  and  the  grid  sizes  were  also  similar 
to  those  used  in  the  base  case  tqxnted  earlier.  Two  cases  were  considered  for  which  where  the  inlet 
Mach  number  of  the  slower,  subsonic  stream  was  held  at  M2  -  0.5.  For  the  supersonic  stream, 
inlet  Mach  number  values  of  M|  »  1.5  and  2.5  were  used.  For  both  streams,  inlet  pressure  and 
density  were  1.01  x  10^  dynes/cm^  and  1.17  x  10"^  gm/cm^  ,  respectively.  The  convective  Mach 
number,  Me  for  the  first  case  is  0.5  and  for  the  second  case  is  1.0. 

Figure  17  shows  the  contours  of  R  within  the  computational  domain  at  a  sequence  of 
times  172  ps  apart  corresponding  to  50  timestq)  intervals  for  the  case  where  Mi  =  1.5.  The  shear 
layer  remains  stable  for  most  of  the  length  of  the  computational  domain  and  then  rolls  up  in  a 
fashion  characteristic  of  substmic  shear  layers .  It  is  of  particular  interest  to  compare  these  results 
to  those  shown  in  Figures  S(b)  and  6.  In  each  case,  M^  =  0.5  and  density  ratios  are  the  same.  The 
velocity  ratios  (U1/U2)  are  1.37  in  Figure  5(b),  1.66  in  Figure  6  and  3.0  for  the  present  case.  The 
computed  values  of  M^^i  and  Mg^  for  a  typical  large  structure  (as  shown  by  the  dashed  line  in 
Figure  16)  is  found  to  be  0.31  and  1.05  respectively.  Thus  even  for  supersonic-subsonic 
interactions,  the  convective  Mach  numbers  for  the  upper  and  lower  stream  are  quite  different  from 
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ttie  theoretical  value. 


In  Figure  18  we  show  the  spectra  of  the  pressure  fluctuations  for  this  case  again  at 
(18.75  cm,  2.0  cm).  Prominent  peaks  are  observed  near  20  kHz  and  in  addition,  there  is  a 
subharmonic  indicating  pairing .  Note  that  for  cases  where  both  streams  are  supersonic,  the  peaks 
are  clustered  near  30  kHz.  The  difference  is  attributed  to  the  existence  of  the  subsonic  stream  and 
acoustic  interactions.  As  the  Mach  number  difference  and  the  convective  Mach  number  are 
increased  ( M}  »  2.5  and  M2  =  0.5) ,  the  pmnt  of  instability  is  found  to  shift  further  downstream. 
The  mixing  between  the  two  streams  also  decreases  within  the  length  of  the  domain  considered.  The 
structures  in  this  case  did  not  exhibit  features  similar  to  those  found  in  subsonic  shear  layers. 
Further  work  (with  improved  outflow  boundary  conditions)  needs  to  be  done  for  characterizing 
subsonic-supersonic  shear  layer. 


V.  CONCLUSIONS 

This  paper  has  presented  numerical  simulations  of  time-dependent  two-dimensional 
unforced  supersoiuc  shear  layers  for  a  splitter-plate  geometry.  The  calculations  were  performed 
using  the  Eulerian  explicit  Flux-Corrected  Transport  algorithm  to  solve  the  convective  transport 
problem.  The  natural  instability  of  supersonic  shear  layer  was  investigated  by  systematically 
varying  die  velocity,  pressure,  and  density  ratios  of  the  two  streams  while  keeping  the  effective  inlet 
momentum  thickness  constant . 

The  suposonic  shear  ayers  were  found  to  be  naturally  unstable  and  to  have  well-defined 
large  structures.  The  organization  and  coherence  of  the  structures  deteriorates  as  the  absolute  Mach 
number  of  the  streams  increases  and  the  Mach  number  difference  between  the  stream  decreases.  For 
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equal  pressure  and  equal  density  streams,  the  mixing  depends  on  the  Mach  number  difference  as 
well  as  the  Mach  number  of  the  faster  stream.  In  this  case,  mixing  decreased  with  decreasing 
convective  Mach  number,  M^;. 

The  most  amplifled  frequency  obtained  by  the  fourier  analysis  of  the  velocity  and 
pressure  fluctuations  depends  on  the  effective  inlet  momentum  thickness  ( the  cell-width  in  the 
cross-stream  direction).  For  the  inviscid  calculations  presented  with  a  step-function  inlet  velocity 
profile,  the  cell  width  in  the  cross-stream  direction  provides  a  length  scale  similar  to  the  initial 
momentum  thickness  of  the  shear  layer.  The  most-amplified  frequency  did  not  depend  on  other 
parameters  like  pressure  and  density  ratios  in  the  shear  layer. 

The  mixing  behavior  of  the  unforced  supersonic  shear  layer  could  not  be  uniquely 
described  by  the  convective  Mach  numba  as  derived  from  the  isentropic  modeP  .  Even  with  the 
same  velocity  and  density  ratios,  different  shear  layers  with  the  same  convective  Mach  number 
show  different  natural  instability  characteristics.  Convective  velocities  M^,!  and  M(;^2 
computed  from  the  predictions  and  were  compared  to  the  (isentropic)  theoretical  value ;  they  were 
found  to  be  very  'different  This  finding  is  in  agreement  with  the  recent  experimental 
observations  by  Papamoschou^.  The  present  results  and  the  experimental  observations^  contradict 
the  current  isentropic  model  of  the  structures  which  predicts  M^j  and  M^^  to  be  equal  or  very 
close. 


Differences  in  pressure  across  a  supersonic  shear  layer  also  enhance  the  convective 
mixing.  When  the  walls  are  close,  shock  waves  originating  at  the  tip  of  the  splitter  plate  reflect 
from  the  nearest  wall  and  intersect  the  shear  layer  before  it  goes  unstable  naturally.  For  • 
equal-temperature  streams  with  the  same  Mach  numbers,  underexpansion  causes  better  mixing  than 
overexpansion. 
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For  equal-pressure  streams,  mixing  is  enhanced  when  the  faster  stream  is  denser 
(cooler).  This  was  found  by  simulating  equal-pressure  shear  layers  with  identical  velocity  ratios  but 
different  density  and  temperature  ratios.  It  is  interesting  to  note  that  Brown  and  Roshko^ 
discovered  the  opposite  trend  for  subsonic  mixing  layers.  This  difference  for  subsonic  and 
supCTsonic  shear  layer  has  been  also  confirmed  by  the  linear  theory 
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Figure  3.  Denssity  contours  from  2.15  x  10’^  gm/cm^  to  2.15  x  10*3  gm/cm^  for  the  base 
cn.se  at  the  end  of  time  step  4000. 
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Figure  4.  Frequency  spectra  in  the  shear  layer  for  the  base  case  at  the  location  18.75  cm  (X), 
2.0  cm  (Y) 
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(c)  Ml  s  4.0,  M2  =  3.5,  Me  =  0.25 


Figure  5. 


Contours  of  mixing  ratio  R  from  0.01  to  0.99  at  the  end  of  lime  step  =  4000  for 
Pi  =  P2  =  1.01  X  lOO  dynes/cm2  ^  pi  =  p2  =  1.17  x  lO'^  gm/cm^. 
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Figure  6.  Contours  of  mixing  ratio  R  from  0.01  to  0.99  at  the  end  of  time  step  =  4000. 

where  Mj  =  2.5,  M2  “  1-5.  =  0.5,  Pj  =  P2  =  1.01  x  10®  dynes/cm^  , 

Pi  =  P2  =  1-17  X  lO'^  gm/cta? 


Figure?.  Contours  of  mixing  ratio  R  from  0.01  to  0.99  at  the  end  of  time  step 
4000  in  a  40  cm  X  3  cm  domain.  Inlet  conditions  are  as  given  in  Figure  2 
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(a)  191  X  61  grid;  time  step  =  1000 


(b)  382  X  61  grid;  time  step  «  1000 


3.0 

g 

0.0 

(c)  382  X  122  grid;  time  step  =  2000 

Figure  9.  Contours  of  mixing  ratio  R  from  0.01  to  0.99  at  a  given  instant  with  different  grid 
densities.  Identical  inlet  velocity  profile  for  each  case. 
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(b)  191  X  61  grid;  time  step  «  4000 


Figure  10.  Contours  of  mixing  ratio  R  from  0.01  to  0.99  at  a  given  instant  with 
different  grid  densities.  Inlet  velocity  profiles  are  different  for  the  two  cases. 
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(a)  Ml  =  5.0,  M2  =  1.5,  Mg  =  0.67,  Pi  =  P2  =1.01  x  10^  dynes/cm2  , 
Pi  s  4.68  X  10‘3  gin/cm^  ,  P2  =  1.17  x  lO"^  gm/cm^ 


X  cm 


(b)  Ml  =  2.5,  M2  =  3.0,  Me  =  0.67,  Pi  =  P2  =  1.01  x  10^  dynes/cm^  , 

Pi  =  1.17  X  10’3  gm/ciTi^  ,  p2  =  4.68  x  10*^  gm/cm^ 

Figure  12.  Contours  of  mixing  ratio  R  from  0.01  to  0.99  at  a  given  instant 
(time  step  =4000)  for  the  same  velocity  ratio  but  different  density  ratios. 
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Ficure  13.  Density  contours  from  1.0  x  lO’^  to  4.0  x  ^ 

underexpanded  case  where  Mi  =  4.0,  M2  =  1.5,  Me  =  1-36,  Pi  =  -0-.  x  10 
dynes/cm*  ,  P2  =  l.pi  x  10®  dynes/cm2  ,  pj  =  2.34  x  lO*-^  gm/cm-^ 

=  1  17  V  in-3  erm/cm3  at  the  end  of  time  step  4000. 
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•  o  n  m  tn  0  99  ai  a  sequence  of  time  steps  (  for  the 
:igu«  14.  Conioars  of  mUing  rano  R  from  0.01  u.  0.99  at  a  s  q 

^  til  ?1011TC  13} 


Figure  15.  Frequency  spectra  in  the  shear  layer  for  pressure  fluctuations  the  case  (shown  in 
Figure  13)  at  the  location  18.75  cm  (X),  2,0  cm  (Y) 
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(a)  Underexpanded  case  where  Mj  =  2.5,  M2  =  1.5,  Me  =  0.633,  Pi  =  2.02  x  10^ 

dynes/cm2,  P2  =  1.01  x  10^  dynes/cm^  and  pj  =  2.34  x  10*3  gm/cm^,  P2  = 
1.17  X  10*3  gni/cm3 


(b)  Overexpanded  case  with  Mi  *  2.5,  M2  =  1.5,  Mg  =  0.369,  Pi  =  2.02  x  10^ 

dynes/cm2,  P2  =  1.01  10^  dyncs/cm2.  and  Pi  =  2.34  x  10*3  gm/cni3,  p2  = 
1.17  X  10*3  gm/cm3 

Figure  16.  Contours  of  mixing  ratio  R  from  0.01  to  0.99  at  a  given  instant 
(time  step  =4(XX))  for  under  and  overexpansion. 
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Frequency  (kHz  x  10) 


Figure  18.  Frequency  spectra  in  the  shear  layer  for  the  pressure  fluctuations  for  the  case  (shown 
in  Figure  16)  at  the  location  18.75  cm  (X),  2.0  cm  (Y) 
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